|
|
|
|
|
|
|
|
EDA Signal Analysis - A Complete Tour |
| Tags | pre-process☁eda☁filter |
Physiological signals are informational resources containing a great potential of knowledge to be extracted, however, this knowledge is not directly accessible after signal acquisition.
In order to convert information into knowledge, physiological signals should be pre-processed, ensuring noise attenuation from unwanted sources and removal of unexpected artifacts. After pre-processing the signal, researcher can be totally secure that his valuable informational resource is ready to be processed, through the extraction of parameters.
Each one of the extracted parameters define an objective way to analyze the acquired and pre-processed physiological data. Subsequent interpretation of the results/extracted parameters will be extremely important to understand which type of events/mechanism are taking place inside the human organism and, at this point, an amazing journey is reaching its fundamental stage with the horizon of achievable discoveries becoming observable by the eyes of knowledge.
| ☌ The current Jupyter Notebook belongs to a set of Notebooks entitled "Complete Tour", where it is presented, in a sequential way, a guide containing recommended steps to be followed, namely: 1) sensor/electrode placement; 2) pre-acquisition notes; 3) pre-processing methods (filtering and motion artifact reduction); 4) detection of specific events and 5) parameter extraction procedures. |
A - Electrodermal Activity (EDA) | Electrode Placement
Taking into consideration the great density of sweat glands (depending on the type of sweat glands reactions such as body temperature control, or emotional/stressful environment feedback can took place) in the palm of the hands and in the fingers
, these become the ideal sites to place the EDA electrodes, as demonstrated in the illustrative example of the following figure:
B - Pre-Acquisition Requirements
In research literature it is considered that the typical informational band of EDA signal is in the range:
-
[0; 3] Hz |
source link
- [0; 10] Hz | source link
-
[0.016; 5] Hz |
source link
-
[0.0167; 0.25] Hz |
source link
-
[0.045; 0.25] Hz |
source link
-
[0.05; 35] Hz |
source link
Accordingly to these articles, informational content of EDA signal reach, at the most broader criteria, the 35 Hz frequency components. So, applying the widely known
Nyquist Theorem
, before starting the data collection, the signal acquisition device must be configured to support sampling rates of at least
70 Hz
, avoiding aliasing problems during the analog-to-digital conversion procedure.
C - Pre-Processing Stage
C0 - Importing Packages and load signal
# biosignalsnotebooks own package.
import biosignalsnotebooks as bsnb
# Scientific programming package.
from numpy import average, array, reshape, sqrt, sort, diff, where, argmax, max
from numbers import Number
# Pacakge dedicated to Wavelet decomposition algorithms.
from pywt import swt, iswt
# Machine-learning dedicated package.
from sklearn.mixture import GaussianMixture
# Gaussian Distribution function.
from scipy.stats import norm
# Built-in packages.
from copy import deepcopy
# Load entire acquisition data.
data, header = bsnb.load("../../signal_samples/eda_hot_surface.txt", get_header=True)
# Store the desired channel (CH3) in the "signal" variable
signal = data["CH3"]
# Sampling rate definition.
sr = header["sampling rate"]
# Raw to uS sample value conversion.
signal_us = bsnb.raw_to_phy("EDA", "bitalino_rev", signal, 10, "uS")
C1 - Filtering Stage 1
In this filtering stage it will be used a conventional bandpass filter, intended to attenuate all signal components outside the typical EDA frequency bands.To compare the influence of multiple filter bandpass frequencies it will be analyzed a restrictive range - [0.045; 0.25] Hz - and a more extended option - [0; 35] Hz.
C1.1 - Generation of filtered signals
# Conventional 1st order Butterworth bandpass filtering
# (due to the restrictive bandwith in order to guarantee system stability it is necessary to reduce the filter order)
# [Option RES]
signal_res = bsnb.bandpass(signal_us, 0.045, 0.25, order=1, fs=sr, use_filtfilt=True)
# Conventional 2nd order Butterworth lowpass filtering
# [Option EXT]
signal_ext = bsnb.lowpass(signal_us, 35, order=2, fs=sr, use_filtfilt=True)
C1.2 - Comparison between filtered and unfiltered signals
With this first filtering stage we achieved a reasonable result:
- The restrictive filtered signal ( signal_res ) was successful in removing movement artifacts, but, unfortunately, it also affect the informational content, as can be seen be the initial concavity and excessively abrupt decrease of the signal after reaching the absolute maximum;
- The extended filtered signal ( signal_ext ) was successful in smoothing the signal but the artifacts and hand movement oscillations were kept;
Taking this information into consideration, it will be necessary to include an additional filtering stage. As a consequence of the preservation of all relevant information, in the subsequent steps only the signal_ext will be analyzed.
The following section is dedicated to minimize the artifacts influence on our EDA signal.
C2 - Filtering Stage 2
After smoothing the EDA signal, removing high-frequency noise components, in the following steps it will be applied an algorithm proposed by Chen et al. ( research article
) to minimize the influence of motion artifacts in the EDA signal.
C2.1 - Decomposition of signal_ext using Stationary Wavelet Transform (SWT)
SWT has the advantage, when comparing with Discrete Wavelet Transform (DWT), of not producing additional "ringing effects in the vicinity of a discontinuity", accordingly to Chen et al. ( research article
).
# SWT 8th level Decomposition using "Haar" mother wavelet (taking into consideration its ability to detect "edges" in the signal) .
swt_orig_coeffs = swt(signal_ext[:32768], "haar", level=8)
# Restriction of filtering algorithm to the coefficients of the last decomposition level.
detail_coeffs = swt_orig_coeffs[0][1]
scaling_coeffs = swt_orig_coeffs[0][0]
C2.2 - Generation of a Gaussian Mixture model
"One Gaussian component describes coefficients centered around zero, and the other describes those spread out at larger values... The Gaussian with smaller variance corresponds to the wavelet coefficients of Skin Conductance Level (SCL) , while the Gaussian with larger variance corresponds to the wavelet coefficients of Skin Conductance Responses (SCRs) " ( research article
)
gaussian_mixt = GaussianMixture(n_components=2, covariance_type="spherical") # "spherical" covariance type ensures that each Gaussian components has its own single variance.
# Reshape data to a column vector format.
detail_coeffs_col = reshape(detail_coeffs, (len(detail_coeffs), 1))
# Fit data to our model object.
gaussian_mixt.fit(detail_coeffs_col)
C2.3 - Determination of the Cumulative Density Function ($\Phi_{mixt}$) of the previously defined Gaussian Mixture
being $weight_1$ and $weight_2$ the relative contribution of each mixture component ($weight_1 + weight_2 = 1$) and $\Phi_1$ and $\Phi_2$ the normal cumulative distribution functions from the original component normal distributions $N(0, \sigma_1^2)$ and $N(0, \sigma_2^2)$.
# Normal distribution function objects.
norm_1 = norm(loc=gaussian_mixt.means_[0][0], scale=sqrt(gaussian_mixt.covariances_[0]))
norm_2 = norm(loc=gaussian_mixt.means_[1][0], scale=sqrt(gaussian_mixt.covariances_[1]))
# Component weights.
weight_1 = gaussian_mixt.weights_[0]
weight_2 = gaussian_mixt.weights_[1]
# CDF values for the coefficients under analysis.
sort_detail_coeffs = sort(detail_coeffs)
norm_1_cdf = norm_1.cdf(sort_detail_coeffs)
norm_2_cdf = norm_2.cdf(sort_detail_coeffs)
# CDF of the Gaussian mixture.
cdf_mixt = weight_1 * norm_1_cdf + weight_2 * norm_2_cdf
C2.4 - Definition of motion artifact removal thresholds using the Cumulative Distribution Function (CDF) of the previously defined Gaussian Mixture model, considering an artifact proportion value $\delta$ equal to 0.01
art_prop = 0.01 # Artifact proportion value.
low_thr = None
high_thr = None
# Check when the CDF mixture function reaches values art_prop / 2 and 1 - art_prop / 2.
for i in range(0, len(norm_1_cdf)):
# Low threshold clause.
if cdf_mixt[i] - cdf_mixt[0] >= art_prop and low_thr == None:
low_thr = sort_detail_coeffs[i]
# High threshold clause.
if cdf_mixt[-1] - cdf_mixt[i] <= art_prop and high_thr == None:
high_thr = sort_detail_coeffs[i]
C2.5 - Removal of wavelet coefficients related with motion artifacts
\begin{equation} d_{2^j}^{2^jk + p} = \begin{cases} d_{2^j}^{2^jk + p}, & \mbox{if } T_{low} < d_{2^j}^{2^jk + p} < T_{high}\\ 0, & otherwise \end{cases} \end{equation}being $d_{2^j}^{2^jk + p}$ the detail coefficient for the $j$ wavelet scale, translation $k$ and shift $p$. $T_{low}$ and $T_{High}$ define the determined lower and higher artifact removal thresholds stored in low_thr and high_thr variables.
filt_detail_coeffs = deepcopy(detail_coeffs)
count_1 = 0
count_2 = 0
for j in range(0, len(filt_detail_coeffs)):
if detail_coeffs[j] <= low_thr or detail_coeffs[j] >= high_thr:
filt_detail_coeffs[j] = 0
else:
continue
# Update of the SWT decomposition tupple.
swt_coeffs = [(array(scaling_coeffs), array(filt_detail_coeffs))]
C2.6 - Signal reconstruction through an inverse SWT
rec_signal = iswt(swt_coeffs, "haar")
C3 - Smoothing stage
In order to remove the constant oscillations in the signal, due to the hand movements, we should execute a smoothing stage through a moving average window. This step needs to be executed only after the removal of the peak artifacts, ensuring that the signal morphology is not disturbed in these points.signal_int = bsnb.smooth(rec_signal, sr * 3) # Moving average witt window size equal to 3 x sampling rate.
# Rescale signal.
signal_int = signal_int * (max(signal_ext) / max(signal_int))
D - Parameter Extraction
In the following cells are presented some extractable parameters that can be extracted from the pre-processed EDA signal, including its formal definition.
| Parameter | Formal Definition | Notes |
| Latency to stimulus onset | $EDR_{lat} = t_{response} - t_{stimulus}$ | $t_{stimulus}$ is the time instant when the stimulus was applied (for example contact of the skin with a hot surface) and $t_{response}$ is the time instant at which the EDA signal starts to change (go out the basal level) |
| EDR amplitude | $EDR_{amp} = EDA_{max} - EDA_{basal}$ | $EDA_{max}$ is the maximum sample value of the EDA signal and $EDA_{basal}$ is the sample value relative to the time instant where the response start ($t_{response}$) |
| EDR rising time (Rise Time) | $EDA_{rist} = t_{max} - t_{response}$ | $t_{max}$ is the time instant of $EDR_{max}$ and $t_{response}$ is the time instant where the response start |
| EDR response peak (Peak Time) | $t_{max}$ | $t_{max}$ is the time instant of $EDR_{max}$ |
| Recovery time to 50% amplitude | $RT_{50\%} = t_{50\%} - t_{max}$ | $RT_{50\%}$ is the time interval that EDA signal takes to decrease 50% of $EDR_{amp}$, $t_{max}$ is the time instant of $EDR_{amp}$ and $t_{50\%}$ is the first time instant, after $EDA_{max}$ where signal amplitude is $EDR_{50\%} = EDA_{max} - 0.50 * EDR_{amp}$ |
| Recovery time to 63% amplitude | $RT_{63\%} = t_{63\%} - t_{max}$ | $RT_{63\%}$ is the time interval that EDA signal takes to decrease 63% of $EDR_{amp}$, $t_{max}$ is the time instant of $EDR_{amp}$ and $t_{63\%}$ is the first time instant, after $EDA_{max}$ where signal amplitude is $EDR_{63\%} = EDA_{max} - 0.63 * EDR_{amp}$ |
# Parameter extraction.
param_dict = {}
#[Latency to stimulus onset]
signal_2nd_der = diff(diff(signal_int)) # Generation of 2nd derivative function.
der_thr = max(signal_2nd_der) # Identification of our major concavity point (start of response time)
response_sample = argmax(signal_2nd_der)
response_time = response_sample / sr
param_dict["Latency to stimulus onset"] = response_time - 0
#[EDR amplitude]
eda_max = max(signal_int)
eda_basal = signal_int[response_sample]
param_dict["EDR amplitude"] = eda_max - eda_basal
#[EDR rising time (Rise Time)]
eda_max_sample = argmax(signal_int)
eda_max_time = eda_max_sample / sr
param_dict["EDR rising time (Rise Time)"] = eda_max_time - response_time
#[EDR response peak (Peak Time)]
param_dict["EDR response peak (Peak Time)"] = eda_max_time
#[Recovery time to 50% amplitude]
time_50 = None
for i in range(eda_max_sample, len(signal_int)):
if signal_int[i] <= eda_max - 0.50 * param_dict["EDR amplitude"]:
time_50 = i / sr
break
param_dict["Recovery time to 50% amplitude"] = time_50 - eda_max_time
#[Recovery time to 63% amplitude]
time_63 = None
for i in range(eda_max_sample, len(signal_int)):
if signal_int[i] <= eda_max - 0.63 * param_dict["EDR amplitude"]:
time_63 = i / sr
break
param_dict["Recovery time to 63% amplitude"] = time_63 - eda_max_time
With the steps described on the current Jupyter Notebook you are able to acquire a broad perspective of the entire methodology involved on EDA signal analysis, starting in the acquisition stage and reaching the parameter extraction.
This is a mere illustrative example aiming to show some possible paths, but, you are free to explore this amazing world by your own.
We hope that you have enjoyed this guide.
biosignalsnotebooks
is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining
Notebooks
!